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SUMMARY 


In this article, we discuss a non-variational V-cycle multigrid algorithm based on the 
cell-centered finite difference scheme for solving a second-order elliptic problem with discontinuous 
coefficients. Due to the poor approximation property of piecewise constant spaces and the 
non-variational nature of our scheme, one step of symmetric linear smoothing in our V-cycle 
multigrid scheme may fail to be a contraction. Again, because of the simple structure of the 
piecewise constant spaces, prolongation and restriction are trivial; we save significant computation 
time with very promising computational results. 

INTRODUCTION 


In the simulation of incompressible fluid flow in porous media, we have to solve at least one 
second-order elliptic equation per each time step. A very important quantity is the Darcy velocity, 
defined by 

u = -KVp U) 

where p is the pressure of the fluid and K is the conductivity. 1C can be written by K. — — , where k 
is a tensor representing the permeability of the medium which can be discontinuous in general, and 
p represents the viscosity of the fluid, p is a continuous function of both time and space variables, 
but may have a very sharp frontal change of values. In other words, p can change rapidly inside the 
interesting domain and the region of rapid change may move as time changes. According to # the 
conservation law of mass balance, the Darcy velocity u must be continuous along the normal 
direction at an element or domain boundary, no matter whether K. is discontinuous or not. 


Now, we consider the 
(p, u) such that 


following simple second-order elliptic equation in mixed form. 

u = -/CVp, in Q = (0, l) 2 C 1R 2 , 

V • u = /, in D, 

p = 0, on dCl, 


Find a pair 


( 2 ) 
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where the conductivity JC(x, y ) = diag(a, b ) is positive and uniformly bounded above and below. 

Because of the discontinuity of /C, the classical solution of p in (2) may not exist. Let (•, •) 
denote L 2 (f2) or (L 2 (fi)) 2 inner product and /f(div; f l) = ju E (L 2 (f2)) 2 | V • u e L 2 (fi) j. We seek 
the solution pair (p, u) E f2) x H(dlv, fi), such that 

(K~ l u,v) = (p, Vv), V v e H{diw,n), 

(V-u, w) = ( f,w ), w E L 2 (fi). 

In [5], error estimates for solving (3) by the cell-centered finite difference scheme are studied, 
with the following results: 

ll-P ~ Pp||l 2 + ||U - 7ru|| L 2 < c/r s ||p|| 1+aifA r, s = 1, 2, (4) 

where V x n is the Raviart-Thomas projection, T are the lines of discontinuity which coincide with 
the grid lines, and (P, U) is the numerical solution of the cell-centered finite difference to 
approximate (3) [5]. Actually, we view the cell-centered finite difference method as a special 
numerical integration of the Raviart-Thomas mixed finite element method [4-6]. For s = 2, (4) is 
the superconvergence error estimate. 

FVom the point of view of mass balance and accuracy, the cell-centered finite difference scheme is 
one of the best numerical schemes to fulfill our goal. In this article, we investigate the efficiency of 
the multigrid algorithm based on the cell-centered finite difference scheme introduced in [5]. 

NUMERICAL SCHEME IN MULTIGRID SETTING 

Let us use the Laplacian operator, —A, to explain the cell-centered finite difference scheme 
stencil. For an interior node, the stencil for —A is (a) in Figure 1. For a corner node, the stencil for 
—A is (b) in Figure 1. For other boundary nodes, the stencil for —A is (c) in Figure 1. For 
discontinuous conductivity, see [5] for details. Now, we consider the uniform grid only. Let Mk 
denote the piecewise constant Raviart-Thomas rectangular pressure space defined on Q with mesh 
size hk = 2“( fc+1 \ k = 0, 1, 2, 3, . . . , J. It is clear that 

Mo C Ml C Mi C .. . C Mj-i c Mj c L 2 (Q). (5) 

With an abuse of notation, for u G Mk, u is either a piecewise function or a vector with its nodal 
values as its entries. On Mk, the cell-centered finite difference approximation is to find P E Mk, 
such that 

A k P = F k = Vkf, & = 0, 1,2, . . ., J. (6) 

Here Vk • L 2 — ► Mk is the L 2 -projection into Mk defined by 

(f,w) = (V k f,w), V tu <E M k , (7) 

and / is the load function of (2). The corresponding stencil of Ak is shown in Figure 1. Our goal is 
to find P E Mj, such that 

AjP = Fj = Vjf. (8) 
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(») (b) (c) 


Figure 1 . Stencils for the Laplacian operator. 


The discrete L 2 -inner product and associated norm on Ai k are denoted by 

{u,v) k = h 2 k v T u and ||u||j| = {u,u) k , u,v £ M k , (9) 

where v T u is the usual algebraic inner product. Let Aj = Aj define an associated bilinear form Aj 
on Mj by 

{Ajw,4>)j = Aj(w, (f>), ( 10 ) 

Before we define A k for 0 < k < J, we first define the prolongation operator I k and the 
restriction operator P k _ v Let I k : M k -\ — * M. k , k = 1, 2, . . . , J be the natural imbedding from 
M. k ~\ to M. k . Thus : M k — * M k -\, the adjoint of I k in (•, •)*, is defined by 

( p k- lW. <f>)k - 1 = (w, Ik4>)k, W £ M k , 4> £ M k -x. (11) 

From (9) and h k _ \ = 2h k , it is clear that in matrix form. Now, we define the bilinear 

form A k _!(-, •) and the matrix A k ^ on M k - 1 for k = J, J - 1, . . . , 2, 1, by 

2A k -i(u, v ) = A k (I k u, I k v), V u, v £ M k -i, (12) 

and the corresponding matrix relation is 

At -i = TAtlt = \Pt.,A k l k . (12') 

Remark 1. It is shown in [5] that for piecewise smooth conductivity tensor 1C, as long as the 
discontinuities coincide with the coarser grid lines 

Ak - i = (l + 0(h 2 )) A k _i. (13) 

In (13) 0(h 2 k ) = Ch\. C depends on the local smoothness of K. but is independent of the jumps. 
Since I k is a simple operator, it is much easier to generate A k _i by (12') than by (6) directly. Of 
course, A k , k = 0, 1, 2, . . . , J — 1, are all positive definite since Aj is, and the spaces are nested. 
Because of (12), our multigrid algorithm can be considered as a black box solver once I k has been 
defined. We mention that (12) holds for three-dimensional problems of —V • (/CVu), with (12') 
being changed to 

At- l = ^A k I i =jpS_ l A i lt. 
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( 14 ) 


We also define the adjoint of I k in A k (-, •), Pk - 1 : Aik —+ AAk-i by 

A k -i(Pk- \u, v) = Ak(% I k v), ueMk , ve M k - 1 - 

To define the smoothing process, we require linear operators R k : M k -* M k for k = 1,2 
These operators may or may not be symmetric with respect to the inner product (•,')*• Let 
A k = D k + L k + Ll, D k be the diagonal part of A ky and L k be the lower triangular part of A k . The 
linear smoothers we have tried are the following relaxation schemes. For 0 < ui < 2, 

(a) Gauss-Seidel: R k = + Lkj an< ^ Pk > 

(b) Jacobi: R k = uD k x , ^ 

Uf 

(c) Richardson: R k = 

where I is the identity operator on M k and A*, is the spectral radius of A k . We allow the relaxation 
parameter u to be different for pre-smoothing and post-smoothing processes in the following 
definition. 

Following [1] the multigrid operator B k : M k - M k is defined by induction and is given as 
follows. The pre-smoother is denoted by R k and the post-smoother by R k . 


V -Cycle Multigrid Algorithm: 

Set B 0 = A^ x . Assume that B k - 1 has been defined and define B k g for g <E M k as follows: 

1. Set x° = 0. 

2. Define x l for l = 1,2,..., m(k) by x e = x e 1 + R k (g — A k x ). 

3. Set y° = + hBk-.Pl, (g - A k x •»«). 

4. Define y e for £ = 1, 2, ... , m{k) by y l = y l ~ l + Rk{g ~ A 

5. Set B k g = y m{k) . 

Remark 2. Since equation (12) holds for all levels, this multigrid algorithm is non-variational 
according to [1], but the approximation property (4) is valid for each level as long as the 
non-variational relation (12) is satisfied. In this algorithm m(k) is a positive integer which may 
vary from level to level. In general this multigrid algorithm is not symmetric in (•, -) k except lor 

Rk = Rl- 

Setting K k = I — R k A k and K k = I - R k A k , it is straightforward to check that 

I-B k A k = KF {k \l ~ hBk-iPl.Ak ]iC (fc) 

= K k 2(k \l ~ hBk-xAk-Sk^Kr^ 
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Equation (16) gives a fundamental recurrence relation for the multigrid operator B k . 


COMPUTATIONAL EXPERIMENTS 


We have tested the multigrid algorithm described in Section 2. We use a power method to 
compute the largest and the smallest eigenvalues of BjAj. 

The linear smoothers we have tried are the following. Let m be a positive integer. m(k) = m for 
all k, 

Si(m) : R k = j-I, R k = R k , 

^(m) : R k — —I, R k = 2 R k , where is the largest eigenvalue of A k , 

Afc 

S 3 (m): R k = (D k + L k )~\ R k = R%, 

S 4 (m): R k = 1.35D k \ R k =±R k , 

Km): Rt =(J± + Lt y\ Rt= ^ + Ll y\ 

S 6 (m): R k =(^D k + L k ') , R k = (2D k + 

Note that only Si(m) and Sz(m) make BjAj Aj ( •, •) symmetric. The rest are neither symmetric 
nor Aj(-, •) symmetric. We also have tried nonlinear smoothers, conjugate gradient, and diagonally 
preconditioned conjugate gradient algorithms. We shall use N(m ) to represent our nonlinear 
multigrid by diagonally preconditioned conjugate gradient smoothers. The reason we choose 
different relaxation numbers comes from the suggestion [3] for an algebraic multigrid algorithm, and 
from our computational experiments. 

We list our test results in Tables 1-9 at the end of this paper for the following problems: 

Ex. 1. Poisson problem: 1C = 1 in (2). 

Ex. 2. Isotropic problem with nearly singular piecewise smooth conductivity: 

K = [o.OOl + 11.1 (l + cos(3.56l7rx) sin(3.00l7rt/))j q , 

_ f 10 -4 , if x > | and y > |, 

^ 1 1, otherwise. 
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Ex. 3. Same kind of problems as Ex. 2: 

K = [O.OOI + 45.l(l + cos(9.431ttx) sin(3.00l7n/))] q, 


9 = 


10 4 , if x > \ and y > |, 
1, otherwise. 


Ex. 4. Anisotropic problem with smooth conductivity: 

K = diag(a, b), 

a = 0.001 + 45.l(l +cos(9.43l7rx)sin(9.43l7ry)), 

b = 0.001 + 45.1 (1 + sin(9.43l7rx) cos(9.43l7ry)) . 

Note that all the solutions of our examples have the superconvergence results proved in [5], i.e., 
satisfying (4) with s = 2. 

In Tables 1 and 6, for example, the second row of Table 1 means J + 1 = 3 level multigrid with 
hj = h A m , Si(l) means A m = min X(BjAj) by Si(l) smoothers, and X M , Si(l) mea.ns 
A M = maxX(BjAj) by S\(l) smoothers. From Table 1, we can see that even when/ Bj j ai s 
to be a reducer, Bj may still be a good preconditioner. In Tables 5-7, it is interesting to see the 
relations of the number of ^-cycles (#V), average contraction numbers (avc) and the time spent on 
the machine (cpu in seconds) when solving a fixed problem on a fixed grid by using different 
multilevels. In Tables 3-5, and 7-9, avc is defined by 

1 A jfaB 

where n = #V is the total number of V^-cycles and ||rjj is the discrete L 2 -norm of the residual 
after the ?th V-cycle. The stop tolerance for all the iterative algorithms is \\r n \\j < e - 10 • Uur 

coarsest grid solver is a diagonal preconditioned conjugate gradient solver with tolerance eo-10 ■ 

In Tables 7-9 “eg” means the standard conjugate gradient algorithm, its corresponding #V 
means the total iteration steps, when \\r n \\ 2 j < « = 10_14 . ^ “ b P c g” means the incom P lete 
factorization preconditioned conjugate gradient algorithm [2]. 
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Table 1. For Ex. 1 



D 

K BBS 

Ajvf , 'S’i(l) 

^1 (2) 

^Af.‘S'i(2) 

4 2 



1.351 

0.788 


8 2 


. 1 

1.804 



16 2 


0.397 

2.394 

0.663 


32 2 


0.367 

3.128 

0.639 


64 2 

5 

0.345 

4.023 

0.623 


128 2 

6 

0.325 

5.106 



256 2 

7 

0.299 

6.417 

0.592 

2.059 


Table 2. For Ex. 1 
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Table 3. For Ex. 1 by 5 j(52(1)) 


Grid 


#F 24 

avc 0.12( 
cpu 35 
#V 26 


J = 2 

J = 3 

34 

45 

0.243 

0.349 

2.1 

1.7 

38 

52 

0.266 

0.384 

11.5 

7.1 

40 

57 

0.27 

0.405 

75.2 

35.2 


Table 4. For Ex. 1 by Bj(Sz(l)) 
J=1|J=21J=31J=4|J=5|J=61J=7|J=8 


Grid J = 1 J — 2 J = 3 J = 4 J = o 
#V 16 22 28 33 34 

64 2 avc 0.050 0.118 0.185 0.256 0.256 

cpu 4.0 1.7 1.5 1-5 1-5 

#V 17 24 31 38 43 


128 2 avc 0.053 0.129 0.204 0.275 0.325 0.345 0.345 

cpu 33.0 8.5 7.0 6.0 6.1 6.4 6.4 

~#y 18 26 34 42 51 58 61 61 

256 2 avc 0.054 0.136 0.219 0.296 0.367 0.417 0.436 0.436 

252.0 61.0 27.0 24.0 28.0 31.0 32.0 




















































Table 8. For Ex. 3 


Grid 


#V 12 
avc 0.028 
cpu 2.3 


#V 14 
avc 0.034 
cpu 9.0 

avc 0.035 
cpu 36.5 


S 5 { 2) I 5 6 (2) S 4 ( 2) bpcg 


11 15 17 26 

0.011 0.047 0.071 

1.0 1.0 1.0 1.2 


17 41 


0.012 I 0.053 0.061 


5 I 3.5 5.5 
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Table 9. For Ex. 4 


Grid 


cpu 32.0 


#V 16 
avc 0.031 
cpu 13.0 


#V 


S 4 (3) 

bpcg 

eg 

21 

27 

313 

0.084 



1.5 

1.3 

2.3 

31 

40 

651 

0.181 



12.0 

5.5 

23.0 

57 

62 

1,329 

0.374 



64.0 

34.0 

161.0 


9.1 x 10 10 


7.2 x 10 u 













































